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We present completely general next-to-leading order predictions for squark and gluino production 
at the LHC, based on the fully automated MadGolem tool. Without any assumptions on the 
mass spectrum we predict production rates and examine the structure of the massless and massive 
quantum corrections. This allows us to quantify theory uncertainties induced by the spectrum 
assumptions commonly made. Going beyond total rates we compare general fixed-order distributions 
to resummed predictions from jet merging. As part of this comprehensive study we present the 
MadGolem treatment of ultraviolet, infrared and on-shell divergences. 
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I. INTRODUCTION 

With the LHC close to completing its 8 TeV run models predicting heavy new particles [I] are under intense 
scrutiny. Experimental searches [2j [3] are probing vast parameter regions in the minimal supersymmetric Standard 
Model (MSSM), most notably those parts of the squark-gluino mass plane which can be described in terms of gravity 
mediation [4 . Inclusive searches for the production and decay of squarks and gluinos plays a leading role and require an 
accurate as well as flexible framework for theory predictions. Understanding decay jets as well as QCD jet radiation [5] 
is a crucial aspect, affecting triggering, rate measurements [5], or kinematic reconstruction. Advanced analysis tools 
like subjet structures [7] increase the need for precise QCD and jets predictions. 

Squark and gluino production to leading order was studied 30 years ago [5]. Next-to-leading order (NLO) QCD 
corrections were first computed almost 20 years ago [SHU] and made public in the Prospino package [T2] 1 . These 
calculations substantially reduce the theoretical uncertainties to the 20 — 30% level. More recently, electroweak 
corrections [15] . resummed predictions [14], and (approximate) NNLO predictions [15] have been made available, 
further decreasing the theoretical uncertainties. Essentially all of these precision studies make simplifying assumptions 
about the squark mass spectrum and focus on improving total cross section predictions. 

In this paper we present numerical results as well as the underlying structure of the completely automized Mad- 
Golem approach to NLO predictions [i"6HT8] . It allows us to go beyond current limitations, like assumptions on the 
supersymmetric mass spectrum or the focus on total rates. We compute the total and differential NLO rates going 
through all squark and gluino pair production channels for several benchmark parameter points. We study in detail 
the structure and numerical impact of the real and virtual QCD and SUSY-QCD effects for each of these channels. 
Particular emphasis we devote to illustrating the reduction of the theoretical uncertainties in total rates and kine- 
matic distributions as a key improvement of NLO predictions. Finally, we conduct a comprehensive comparison of the 
fixed-order differential cross sections with those obtained by multi-jet matrix element merging, including a variation 
of the renormalization and factorization scales. Many details on the computation and its numerical validation are 
included in the four appendices. 

This study, alongside with its earlier squark-neutralino 16J and sgluon-pair |17j counter parts, are examples of fully 
automized NLO computations in TeV-scale new physics models. This kind of automation will significantly enhance 
the availability of precision predictions for LHC observables in and beyond the Standard Model [T5], at a time when 
standard new physics scenarios for the LHC are becoming less and less likely. MADGOLEMis an independent, highly 
modular add-on to MadGraph [20j[2T], benefiting from its event simulation and analysis features. It generates all 
tree-level diagrams and helicity amplitudes in the MadGraph v4.5 framework [2U] and relies on the Helas [2"2"] 
library for the numerical evaluation. The one- loop amplitudes are generated by Qgraf [33J and Golem [2~4l I25j. 
Supersymmetric counterterms and Catani-Seymour dipoles [26] are part of our model implementation and can easily 
be adapted for other new physics models. The subtraction of infrared and on-shell [HI [27] divergences is completely 
automized. MadGolem is currently undergoing final tests and will be released to the LHC community soon. 

II. RATES 

As a starting point we systematically analyze squark and gluino production rates at next-to-leading order. We 
entertain all possible production channels at the LHC involving pairs of squarks and gluinos in the final state: 

pp^ qq,qq*,qg,gg . (1) 

Following the typical decay signature we focus on the dominant first and second generation squarks q = 
ul.r, Al,r, sl,_r, cl._r- The associated quarks we can safely assume to be massless. Moreover, we disregard flavor- 
mixing, i.e. the SUSY-QCD couplings are flavor-diagonal. Further removing this latter assumption is foreseen in 
the MadGolem setup. In our numerical analysis we use the CTEQ6L1 and CTEQ6M parton densities [28] ■ Un- 
less stated otherwise, we fix both the central renormalization and factorization scales to the average final-state mass 
Mi? = Mf = M G = {mi+m2) /2. From previous studies, this choice is known to lead to perturbatively stable results [11] . 

As real corrections we include all channels with a three-particle final state, in which a light parton accompanies the 
heavy superpartners. The associated infrared divergences we subtract using Catani-Seymour dipoles [26] . generalized 



All results shown in this paper have been checked to agree with PROSPIN02 wherever applicable. 
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to include the massive colored SUSY particles. Details on this implementation are included in Appendix [XJ With the 
help of an FKS-like cutoff a [29] we can select the phase space regions covered by the dipole subtraction to include 
more (a — > 1) or less (a <§; 1) of the non-divergent phase space. Our default choice is a = 1, but the total rates must 
not change with varying a. 

Virtual corrections include the one-loop exchange of virtual gluons and gluinos. The standard 't Hooft-Feynman 
gauge is employed for internal gluons to avoid higher rank loop integrals. Accordingly, Fadeev-Popov ghosts appear in 
the gluon self-energy and in the three-gluon vertex corrections. Ultraviolet divergences are cancelled by renormalizing 
the strong coupling constant and all masses. Supersymmetry identifies the strong gauge coupling constant g s and the 
Yukawa coupling of the quark-squark-gluino interaction, g s . At the one- loop level dimensional regularization induces 
an explicit breaking of this symmetry via the mismatch between the 2 fermionic gluino components and the (2 — 2e) 
degrees of freedom of the transverse gluon field. We restore the underlying supersymmetry with an appropriate finite 
counter term [111 130] . Details on the renormalization procedure can be found in Appendix [Dj 

Finally, we have to remove potential divergences in the three-body phase space due to intermediate resonant 
states. An example is the appearance of on-shell gluinos as part of the correction to squark-antisquark production, 
qg — > qg — > qq*q |11U27| . In addition to the technical complication of a divergent rate these on-shell states introduce a 
double counting once we sum all squark and gluino production rates to next-to-leading order. In the Standard Model 
a similar problem appears in Wt single top production which requires a separation from top pair production, so our 
MadGolem implementation should benefit Standard Model processes as well. 

Following the Prospino scheme we remove on-shell divergences locally through a point-by-point subtraction over 
the entire phase space. Off-shell pieces in the limit of vanishing particle width are genuine parts of the NLO real 
emission and hence left untouched. This procedure preserves the gauge invariance of the entire matrix element as 
well as the spin correlations between the intermediate particles and the final state. The subtraction terms have a 
Breit-Wigner shape and are automatically generated. 

Note that for an actual observable we of course need to combine the pair production and associated production 
channels. Initial-state jet radiation at the LHC may be as hard as decay jets [5] and thus cannot be distinguished on 
an event-by-event basis. A detailed account of our on-shell subtraction is provided in Appendix |B[ 



A. Parameter space 



The effect of NLO corrections on LHC cross sections varies from production channel to production channel and from 
one mass spectrum to another. In particular the hierarchy between squark and gluino masses affects the behavior 
of QCD corrections. In Tab. [i] we list a set of conventional mass spectra [31] which in the following we will use 
to study squark and gluino pair production in some detail. Scenarios labelled CMSSM- # (constrained MSSM) are 
derived from GUT-scale universality conditions with squark and gluino masses above 1 TeV. Each of them exhibits 
a different squark-gluino mass hierarchy. The benchmark points denoted as mGMSB-# and mAMSB-# represent 
gauge-mediated and anomaly-mediated supersymmetry breaking. While the entire spectrum including the weak 
gauginos is very different compared to CMSSM-type scenarios, they are considerably less distinctive when we limit 
ourselves to light-flavor squark and gluino masses. The position of each of the benchmark points in the squark- 
gluino mass plane defines a maximal mediation scale for SUSY breaking, but this argument cannot simply be turned 
around [I]. 



In Tabs. [II] and III we document a comprehensive numerical survey over these MSSM parameter for LHC center- 
of-mass energies of ^/ S = 14 TeV and >/S = 8 TeV. These cross sections can be used to test MadGolem once it is 
publicly available. Using the general MadGolem setup it is possible to separate the squark flavor and chirality in 
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Table I: Squark and gluino masses (in GeV) for different benchmark points. 
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CMSSM 10.2.2 
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3.0 
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2.5 
27.5 
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1.5 



5.8 
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29.2 
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1.44 
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4.6 
3.7 
1.9 
21.1 
17.1 
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6.0 
4.9 
2.6 
27.8 
22.5 
3.0 



1.30 
1.32 
1.33 
1.32 
1.32 
1.32 



16.0 
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7.7 
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CMSSM 10.2.2 
CMSSM 40.2.2 
CMSSM 40.3.2 
mGMSB 1.2 
mGMSB 2.1.2 
mAMSB 1.3 



78.7 
93.5 
159.4 
467.0 
777.0 
54.4 



108.6 
131.3 
239.5 
610.6 
1077.6 
78.1 



1.38 
1.40 
1.50 
1.31 
1.39 
1.44 



87.7 120.3 

101.7 142.3 

165.6 248.2 

511.4 665.4 

868.0 1193.9 

53.5 77.0 



1.37 
1.40 
1.50 
1.30 
1.38 
1.44 



2.3 
2.8 
5.2 
18.7 
33.6 
1.5 



3.8 
4.6 
9.0 
28.3 
52.5 
2.6 



1.63 58.2 

1.65 68.7 

1.73 116.3 

1.52 371.2 

1.56 638.1 

1.71 36.3 



83.6 
100.5 
182.0 
503.3 
914.6 

54.5 



1.44 
1.46 
1.57 
1.36 
1.43 
1.50 



23.3 53.4 2.29 

41.1 94.5 2.30 

249.2 552.9 2.22 

222.8 453.4 2.03 

849.6 1755.0 2.07 

19.0 46.1 2.42 



Table II: Total cross sections (in fb) and corresponding K factors for squark and gluino production at \/~S = 14 TeV. The 
renormalization and factorization scales are chosen as the average final state mass. The notation ud indicates the summation 
over all possible final-state chiralities ud = Uhdh + u^dR + ur<1l + ur<1r. Symmetry factors 1/2 are automatically included, 
when applicable. 
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Table III: Total cross sections (in fb) and corresponding K factors for squark and gluino production at \fS = 8 TeV. All scales 
are chosen at the average final state mass. 



squark pair production and in associated squark-gluino production. The size of the NLO QCD effects we express 
through the consistent ratio K = ct nlo /(t lo , in spite of some well-known problems with the convergence of the LO 
parton densities we will notice below. 

Already at the level of total cross sections we confirm a number of well known general trends which are essentially 
common to all production mechanisms. First, the significance of the QCD quantum effects manifests itself as sizable K 
factors spanning the range K <~ 1.1 — 2.4 for 14 TeV center-of-mass energy. For 8 TeV we observe uncomfortably large 
K factors which have nothing to do with real or virtual QCD corrections. Instead, they indicate poor perturbative 
behavior of the CTEQ parton densities. These suppress LO production rates of heavy particles with O(TeV) masses, 
while the NLO predictions are perturbatively stable. This can be checked by comparing the CTEQ rate predictions 
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to other parton densities. 

Second, the different color charges of squarks and gluinos as well as their spin are clearly reflected in the production 
rates. Interactions among color octets will give larger rates than color triplets. Similarly, fermion pairs yield larger 
cross sections than scalar pairs. This effect is not only observed in the LO and NLO rates but also in the relative K 
factor. 

Third, both the total LO and NLO cross sections decrease with increasing superpartner masses. Non-trivial effects 
can for example be understood from the threshold behavior of virtual corrections and the real emission, which may 
in part overcome the phase space suppression of the NLO diagrams [TT] . We will expand on all these aspects later in 
this Section. 



B. Squark pair production 

Squark pair production can lead to a multitude of final states, which we first classify into two basic categories: 

1. squark-squark pairs qq, to leading order mediated by ^-channel gluino interchange between colliding quarks. This 
mechanism is flavor-locked, so first generation squarks will dominate. In particular in proton-proton collisions at 
large parton- a: values this channel will contribute large cross sections because it links incoming valence quarks. 

2. squark-antisquark pairs qq* with three distinct sub-channels: qq annihilation through an s-channel gluon; qq 
scattering via a f-channel gluino, and gg fusion with s-channel and ^-channel diagrams. Due to the large adjoint 
color charge and the higher spin representations involved the gg initial-state dominates at the LHC up to 
moderate parton-x values. In the absence of flavor mixing, the gluino-induced sub-channel is flavor-locked to 
the initial state while the other two are flavor-locked within the final state. First and second generation squarks 
will therefore contribute with similar rates. All but the gluino-induced channels will also lead to sbottom and 
stop pair production [9]. 

The predicted LO and NLO rates alongside their K factors we document in Tables [XT] and |III| The production of 
squark pairs qq yields cross sections of 10 to 100 fb for squark and gluino masses around 1 TeV. The squark-antisquark 
rates for this mass range are roughly one order of magnitude smaller. These cross sections are highly sensitive to 
the strongly interacting superpartner masses. This is largely due to kinematics, i.e. the different squark masses in 
each benchmark point. For instance, the production of the lighter right-handed squarks comes with larger production 
rates than that of their left-handed counterparts. According to Tab. [I] this is true for all benchmark points except for 
mAMSB 1.3. This means that in a squark-(anti)squark sample right-handed squarks will be overrepresented. This 
can be a problem if the NLO computation does not keep track of the different masses of left-handed and right-handed 
quarks. 

In contrast, we see that the K factors barely change between benchmark points, because the bulk of the NLO 
effects are genuine QCD effects. However, all K factors range around K ~ 1.2 for squark-squark production - 
correspondingly, for squark-antisquark production they render K ~ 1.2 — 1.5 depending on the specific channel. 
This effect is used by Prospino2.1, where the different squark masses only enter at leading order, while the NLO 
corrections are computed with a universal squark mass. Some sample Feynman diagrams we show in Fig. [T] The 
supersymmetric QCD corrections including one-loop squark and gluino loops are power suppressed by the heavy 
particle masses. 

An interesting observation we make for squark pairs with different chiralities, e.g. u^ur. As mentioned above, 
all qq channels proceed via a i-channel gluino. For identical final-state chiralities, the gluino propagator corresponds 
to a mass insertion - enhancing the LO rates for heavy gluinos. This is not true for ulUr production, where we 
probe the fl term in the gluino propagator. This difference can be read off Tab. ITT| The ul^l & n d urUr channels are 
suppressed from CMSSM 10.2.2 to CMSSM 40.3.2, following a decrease in the gluino mass. On the other hand, the 




Figure 1: Sample Feynman diagrams for squark-antisquark production to NLO. Virtual corrections involve the exchange of 
gluons, gluinos and squarks. Real corrections denote the emission of one quark or gluon. 




Figure 2: Cross sections for ulu* l production for the different initial states as a function of the squark and the gluino masses. 
The qq process (left) includes also the qg crossed-channels. Together with mi 1 we vary all squark and gluino masses such that 
the mass splittings of the CMSSM 10.2.2 benchmark point are kept. In the lower panels we evaluate the relative size of the 
NLO cross section with respect to the total LO rate for each sub-channel. 



ulUr rate remains quite constant. 
Kll ~ K RR < K LR . 



This different behavior is also visible from their K factors, which are ordered as 



In Fig. [2] we separate the real and virtual QCD and SUSY-QCD corrections for ulu* l production as a function of 



the final state mass 



All the other heavy masses we vary simultaneously, keeping the absolute mass splittings of 



the CMSSM 10.2.2 benchmark point shown in Tab.|T] The two main partonic subprocesses contributing to the process 
we show separately. The separation into real and virtual corrections we define through Catani-Seymour dipoles with 
a FKS-like cutoff a — 1. The integrated dipoles count towards the virtual corrections while only the hard gluon 
radiation counts towards the real corrections. This is the reason why the real corrections appear negligible. The cross 
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Figure 3: Relative shift Act nlo /ct lo for the different parts of the virtual corrections to qq/gg — > UtitJ, production. All squark 
and gluino masses we vary in parallel, just like in Fig. [2] 
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MjirV H F V V R /V° M*V 

Figure 4: Renormalization and factorization scale dependence for squark pair production pp — > ulul (upper) and pp — > ulu* l 
(lower). The plots trace a contour in the [j,r-[j.f plane in the range /i = (0.1 — 10) x fi° with fi = m% L - All MSSM parameters 
follow the CMSSM 10.2.2 benchmark point in Tab. |T] 



sections for both the gluon fusion gg and the quark-antiquark qq subprocesses are essentially determined by the squark 
masses and the corresponding phase space suppression. The gluon fusion dominates in the lower squark mass range, 
contributing with rates of roughly a factor 2 above the qq mechanism. Conversely, the gg channel depletes slightly 
faster than the qq, especially for large squark masses. This can be traced back to the respective scaling behavior 
of the cross sections [S] as a function of the partonic energy, and its correlation to the parton luminosities. Indeed, 
heavier final-states probe larger parton- x values — this being the region where the quark parton densities become 
more competitive, while the gluon luminosity depletes. 

The lower panels of Fig. [2] show the relative size of the NLO contributions with respect to the total LO rate. While 
cr vlrtual /er LO grows with increasing squark masses, specially for the gg initial state, cr rcal /<7 LO stays constant. This 
effect is related to threshold enhancements: first, a long-range gluon exchange between slowly moving squarks in the 
gg — > uu* channel gives rise to a Coulomb singularity a ~ na s /f3, where /3 denotes the relative squark velocity in 

the center-of-mass frame, (3 = \J 1 — Am\/S. This is nothing but the well-known Sommcrfeld enhancement 32J. The 
associated threshold singularity cancels the leading a ~ (3 dependence from the phase space and leads to finite rates 
but divergent K factors [TT]. In addition, there exists a logarithmic enhancement a ~ [Alog 2 (/3) -I- i?log(8/3 2 )] from 
initial-state soft gluon radiation. This second effect is common to the gg and qq initial states. Threshold effects can 
be re-summed to improve the precision of the cross section prediction [14] . 

The internal architecture of the virtual corrections we analyze in Fig. [3j Virtual diagrams come in different one- loop 
topologies: self-energy and wave-function corrections, three-point vertex corrections, and box corrections. The box 
diagrams also include the one-loop corrections to the quartic ggqq vertex. Again, we assume the specific flavor/chirality 
final state ulu* l with the CMSSM 10.2.2 parameter point. Just like in Fig. [2] the masses vary in parallel, keeping 
the splitting constant. The threshold effects discussed in the previous paragraph are nicely visible in the increasing 
ratio A(7 NLO /<7 LO for the boxes and the integrated dipoles, where the quantity Act nlo /<t nlo accounts for the genuine 
0(a s ) NLO contributions. This enhancement leads to sizable quantum effects in the 30% — 70% range for the gg 
initial state. 

For the qq-initiated subprocess the integrated dipoles are numerically far smaller. The bulk of the virtual corrections 
is driven by the boxes, the gluino self-energy, and the negative quark-squark-gluino vertex correction. Their remark- 
able size we can trace back to mass insertions in the gluino-mediated diagrams. Barring these dominant sources, 
Fig. [3] illustrates that all remaining NLO contributions stay at the ~ 5% level or below. In the absence of threshold 




Figure 5: Cross sections a(pp — > ulul) (left) and u(pp — > ulu*l) (right) as a function of the squark mass. The band 
corresponds to the scale variation envelope /i°/2 < (ir,f < 2fjP, where /x° = rriu L . The central MSSM parameters are given by 
the CMSSM 10.2.2 benchmark point. The squark and gluino masses we vary in parallel, just like in Fig. [2] 



effects, all these pieces are insensitive to the squark mass. As a consequence, both the LO and the NLO cross sections 
undergo essentially the same phase space suppression as a function of the final state mass. Because we vary all masses 
in parallel this is also indicative of the dominance of the gluon-mediated QCD effects as compared to SUSY-QCD 
corrections. In the large-mass regime the latter have to be power suppressed, matching on to the decoupling regime. 

The fact that cross section predictions increase, i.e. exclusion limits become stronger once we include NLO cross 
sections is only a superficial effect of the improved QCD predictions. The main reason for higher order calculations 
is the increased precision, reflected in the renormalization and factorization scale dependence. As is well know, these 
scale dependences do not have to be an accurate measure of the theoretical uncertainty. This can be seen for example 
in Drell-Yan-type processes at the LHC where the LO factorization scale dependence hugely undershoots the known 
NLO corrections. For the pair production of heavy states mediated by the strong interaction the detailed studies of 
top pairs give us hope that the scale dependence can be used as a reasonable error estimate. 

In Fig. [4] we trace the scale dependences of squark-squark and squark-antisquark production. Note that such a 
separate scale variation is not possible in Prospino, where both scales are identified in the analytic expressions. We 
profile the behavior of a (/z) and a N (/z) for an independent variation of the renormalization and the factorization 
scales in the range /i°/10 < ^b.,f < IO/1 . As usual, the central scale choice is /jP = rriu L . The path across the -^f 
plane we illustrate in the little square in the left panel. The numerical results are again given for the CMSSM 10.2.2 
parameter point and \/S = 14 TeV. As expected, the renormalization scale dependence dominates the leading order 
scale dependence. Unlike in other cases there is no cancellation between the renormalization and the factorization 
scale dependences. The stabilization of the scale dependence manifests itself as smoother NLO slope. While the LO 
scale variation covers an 0(100%) band, the improved NLO uncertainty is limited to O(30%). Interestingly, the NLO 
plateau at small scales is not generated by a combination of the two scale dependences, but is visible for a variation 
of the renormalization scale alone at fixed small values of the factorization scale. 



In Fig. [5] we show the usual LO and NLO cross sections as a function of the final-state mass m^ L . The error bar 
around the central values represents a simultaneous scale variation [/i°/2, 2/i ]. Both error bands nicely overlap and 
reflect, for ulul, a reduction of the theoretical uncertainties from O(50%) at LO down to O(20%) at NLO - similarly, 
from 0(60%) down to 0(30%) for u L u* L . 



The most significant upgrade of the MadGolem automated framework compared to previous calculations is that 
we do not have to assume any simplifying relations between the supersymmetric masses. We can freely sweep over the 
entire parameter space of a given model, varying each input parameters independently. This differs from Prospino or 
other precision tools which rely on a single mass scale for all light-flavor squarks for all next-to-leading order effects. 
A fully general scan as shown in Tabs. TlpTT is thus beyond the reach of these tools. 



Figure [6] shows quantitative results for this generalized NLO computation. As an example we focus on the (partially 
inclusive) production of all first-generation squark pairs pp — > qq, with q — ul, ur, d,L,d,R and examine an independent 
variation of the different squark masses. In the left two panels we study the effect of a right-left mass separation while 
identifying s-up and s-charm masses as well as s-down and s-strange masses to the CMSSM 10.2.2 and mGMSB 2.1.2 
values shown in Tab. [TJ We show the change in the total squark pair cross sections with a growing mass splitting 
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Figure 6: Cross sections for squark pair production pp — » qq (q = UL,R,dL,n) as a function of mass splittings. In the left 
panels we vary the right-left splitting keeping the flavor splitting constant. In the right panels we vary the u-d flavor splitting 
fixing the right-left splitting. We show the shift with respect to the degenerate spectrum with the masses and the total rates 
Co = cr(Am = 0) given in each panel. 



Amj;„j, = mq R — rriq L , where q = u,d, c, s. The results we evaluate in terms of \a — Co \/cr , where Co denotes the cross 
section for Aiur-l = 0. The associated K factor is displayed at the bottom of each panel. In the right panels we 
show the same analysis for an up-versus-down squark mass splitting Am u -d = ma ~ m ^ with u — u,c and d — d,s. 
The masses of their respective chiral components are separated as in the corresponding benchmark scenarios. 

The results highlight, first of all, that a fully flexible mass spectrum leaves a measurable footprint in the total 
cross sections. The rates change by 0(5 — 20%) for a squark mass splitting of 10 — 100 GeV, as commonly featured 
by MSSM benchmark points. These effects lie roughly in the same ball park as higher order corrections beyond the 
fixed next-to-leading order predictions. At the same time, the K factors stay essentially constant with a varying mass 
splitting. This follows from the fact that while kinematic effects change the cross sections significantly, the NLO 
corrections are mostly sensitive to mass splittings through SUSY-QCD effects which are typically mass suppressed. 
For squark pair production, in particular, the NLO rate dependence on Attir-l and on Am u -d is affected by the 
gluino self-energy corrections. Our MadGolem results confirm that taking into account the full mass spectrum in the 
LO rate predictions and re-weighting them with a K factor computed most efficiently for a mass degenerate spectrum 
gives an accurate estimate of the full NLO rates. 




Figure 7: Cross sections for a(pp — > ul§) as a function of the squark mass mu L - The band corresponds to a scale variation 
H°/2 < flR 3 F < 2/x°, where fi° = (ma L + mg)/2. The MSSM parameters are given by the CMSSM 10.2.2 benchmark point. 
The squark and gluino masses we vary in parallel, just like in Fig. [2] 
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Figure 8: Renormalization and factorization scale dependence for ul§ associated production. The plot traces a contour in the 
(AR-fJ-F plane in the range fi — (0.1 — 10) x fi° with fj? = (mo L + nig)/2. All parameters are the same as for Fig. [4] with mass 
values rriu L = 1162 GeV and m§ = 1255 GeV. 




Figure 9: Cross sections for squark-gluino production pp — > qg (q = UL,R,dL,R,) as a function of mass splittings, within the 
same setup as in Fig. [6] 

C. Squark-gluino production 

Unlike squark or gluino pair production the associated production process does not have a QCD-only component 
and is always flavor locked, qg — > qg. This makes it the most model dependent signature. First generation squarks, 
mostly UL t R, will be copiously produced, and some of the structures will be reminiscent of the electroweak production 
of squarks with electroweak gauginos, pp —> q\ [To] - 

Moreover, in this particular channel on-shell divergences can have a twofold origin: they can either stem from 
an on-shell gluino or an on-shell squark, depending on which of these particles is heavier. This makes associated 
squark-gluino production the key channel to test our numerical MadGolem implementation of automized on-shell 
subtraction. 

Several qualitative expectations we can nicely confirm from Tabs.[IT|and |IH| For instance, we see how u^g production 
dominates over the charge conjugated channel u* L g, simply due to the valence u quark. This is also the reason why 
the QCD corrections are larger for the u* L g process, because gg-initiated NLO contributions are not suppressed by 
the relative size of the underlying parton luminosities. 

The dependence on the final state masses we show in Fig. [7j where we display the total cross sections a(pp — > u^g) 
as a function of the final-state squark mass ma L , noting that the gluino mass is changed together with the squark 
mass. The total cross section is pulled down by roughly three orders of magnitude when the final-state mass increases 
by a factor of three. We find cross sections as large as <t(ul<?) ~ O(10) pb for m Si < 500 GeV, which fall down 
to O(10) fb for m u ~ L < 1.5 TeV. A remarkable reduction in scale dependence of 0(60)% down to 0(20)% can be 
assessed by contrasting the LO and NLO uncertainty bands. A complementary viewpoint we provide in Fig. [HJ where 
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Figure 10: Renormalization and factorization scale dependence for gluino pair production. The plot traces a contour in the 
fJ-R-fJ-F plane in the range \i — (0.1 — 10) x fj? with fj? = m§. All parameters are the same as for Fig. HI with m§ = 1255 GeV. 



we probe scale variations of the total cross section as usually in the two-dimensional renormalization vs factorization 
scale plane. 

Finally, we address the effect of a general squark mass pattern on the total rates. Our analysis follows Pig.[6j now for 
the process pp — > qg with q = ul,r, cIl,r- Again, we study the different CMSSM 10.2.2 and mGMSB 2.1.2 scenarios. 
For each of them, we explore the relative change in the total rate \a — a°\/a° when we increase mass splittings from 
zero ((Tq)- We separately examine i) fixing all left-handed and right-handed squarks at one common mass value and 
increasing the right-left mass splitting Atur-l', and ii) setting a common mass for up- type and down- type squarks 
and increasing Am u _ d . Similarly to the squark pair case, in Fig. [9] we find variations up to 20% for mass splittings 
up to 0(100) GeV. The LO and NLO cross sections scale in parallel, with minor differences at the per-cent level. 
As expected from the squark pair case the footprint of a non-degenerate squark spectrum factorizes from the QCD 
corrections, with remaining non-factorizing effects the level of a few per-cent. 



D. Gluino pair production 



Finally, similar phenomenplogical trends we identify for gluino-pair final-states. The NLO effects are particularly 
sizable (cf. Tables O and III I with K factors in the ball-park of ~ 2 for yfS = 14 TeV, and even surpassing K ~ 2.5 
for the lower nominal LHC energy yfS = 8 TeV. These results essentially reproduce what is included in Prospino. 

The separate dependence on the factorization and renormalization scales we display in Fig. 10 As before, the 
simultaneous scale variation captures the complete theoretical uncertainty well. Including the NLO corrections sig- 
nificantly reduces the dependence on the renormalization as well as on the factorization scale. In Fig. [TT] we show 
the envelope of the scale variation together with the central LO and NLO rate predictions as a function of the gluino 




Figure 11: Cross sections for o(jpp — > gg) as a function of the gluino mass mg. The band corresponds to a scale variation 
fX°/2 < g.R.F < 2fi° with fj? = rrig. The MSSM parameters are given by the CMSSM 10.2.2 benchmark point. The squark and 
gluino masses we vary in parallel, just like in Fig. [2] 
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Figure 12: Normalized transverse momentum distributions for different processes for the benchmark points CMSSM 10.2.2 and 
mGMSB 2.1.2. We compare NLO predictions to LO jet merging [33] with three different setups: up to one hard gluon; up to 
two hard gluons; up to one hard quark or gluon jet. The latter two we only display when differences are visible. 



mass. It reflects a reduction of the theoretical uncertainties from 0(70%) at LO down to 0(30%) at NLO. In spite of 
the large K factor triggered by the LO parton densities the two bands nicely overlap. 

From the discussions of the squark-pair and squark-gluino channels we expect the effect of the non-degenerate 
squark spectrum on gluino pair production to be small. The leading phase space effects which appear as factorizing 
corrections in the other production channels are absent for gluino pairs. Only the LO squark mass dependence through 
the f-channel exchange diagram and the specific NLO loop effects remain. For the mGMSB 2.1.2 benchmark point 
we compute these effects and find a deviation of roughly 2% when we allow for a right-left mass splitting of 100 GeV. 



III. DISTRIBUTIONS 



All through Sec. [IT] we have limited our discussions to total cross sections. This corresponds to the way higher-order 
corrections to new physics processes are usually implemented in experimental analyses. Event simulation including 
all differential cross section is performed by any of the parton-shower Monte Carlos. Because the hard process scale 
is given by the heavy particle masses it is usually well above the typical jet momenta required by inclusive searches. 
This means that the parton shower approximation is justified [5 while the total cross sections have to be corrected 
for higher-order effects. In the original Prospino calculations transverse momentum and rapidity distributions 
for the heavy squarks and gluinos were studied [9HTT]. indicating that no large NLO effects should be expected. 
MadGolem allows us to include a comprehensive study of distributions in this paper. 



A. Fixed order vs jet merging 



To make quantitative statements beyond total cross sections we use MadGolem to compute NLO distributions for 
different squark and gluino final states. Because the MadGolem output is weighted events for the regularized virtual 
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Figure 13: Two-dimensional distributions for squark pair production pp — > ulUl at "v/S* = 14 TeV as contour plots in the 
Pt(ul)-u(ul) plane. The different panels show the results from LO (left), NLO (center), and jet merging (right). While the 
LO result is shown to scale the two right histograms are normalized to unity. We use the CMSSM 10.2.2 parameters. 



and real emission channels we can plot any distribution which makes sense in perturbative QCD. The only limitation 
is the validity of fixed-order QCD, reaching its limitations for example when studying the jet recoil against the heavy 
squark-gluino system. 

For comparison we do not rely on the usual parton shower simulations, but employ more modern matrix element 
and parton shower merging |33H35j . We generate tree-level matrix element events with zero, one, or two hard jets 
with the help of MadGraph5 [21 and combine them with each other and with the Pythia [35] shower using 
the MLM procedure [33J as implemented in MadGraph. When defining the hard matrix element corrections we 
follow three different approaches. First, we include up to one additional hard gluon in the matrix elements. This 
automatically excludes all topologies which could lead to on-shell divergences. Second, we instead allow for two 
additional hard gluons in the matrix elements. As before, we avoid any issues with on-shell singularities. Finally, we 
generate samples with one additional quark or gluon. In this case, on-shell divergences will appear just like for the real 
emission contributing to the NLO rate. These singularities we remove using the numerical prescription implemented 
in MadGraph [3T]. It subtracts all events with phase space configurations close to the on-shell poles. While this 
subtraction is not equivalent to the consistent Prospino scheme and does not have a well-defined zero-width limit 
we have checked that it gives numerically similar results as long as we only compare normalized distributions. 

Our results for the transverse momentum distributions of squarks and gluinos we present in Fig. |12| We focus on 
the CMSSM 10.2.2 and mGMSB 2.1.2 benchmark points as representative MSSM scenarios. They exemplify both 
possible squark-gluino mass hierarchies. As described above, we show the NLO predictions and the one-gluon merged 
results. Comparing different jet merging setups we confirm that adding a second hard gluon does not change our 
results beyond numerical precision, so we do not show it separately. This is an effect of the large hard scale in the 
process which renders the parton shower for the second radiated gluon an excellent approximation. Results allowing 
for one additional quark or gluon jet we only show when the curves are visibly different from the one-gluon case. 
The general agreement of all three merging results shows that once the double counting from the on-shell states is 
removed the bulk of the NLO real emission comes from gluons. Moreover, this gluon radiation is well described by 
the Pythia parton shower, as long as the produced particles are heavy [5], However, once an experimental analysis 
becomes particularly sensitive to the jet recoil it might pay off to check the parton shower results with a merged 
sample 0137]. 

The comparison with the NLO prediction shows that the usual assumption about the stability of the main dis- 
tributions is indeed correct. The normalized distributions from the fixed-order NLO calculation and from multi-jet 
merging agree very well. As alluded to above, the multi-jet merging predictions in turn agree well with the parton 
shower. In spite of the remarkable agreement between both descriptions, some mild departures arc visible. We can 
essentially understand them as a fingerprint of the extra recoil jets involved in the matched samples. For example, in 
some cases the jet-merging predictions become slightly harder than the NLO results because they take into account a 
second radiated jet. On the other hand, the squarks and gluinos we are studying are so heavy that it is unlikely that 
jet radiation makes a big difference to them. 



is corn- 



In MadGolem the generation of any kind of fixed-order distributions, such as those displayed in Fig. 12 
pletely automated. This constitutes a substantial improvement for precision BSM predictions. Distributions can be 
computed for a single kinematic variable, but also two-dimensionally. For example, we show the NLO phase space 



dependence on the transverse momentum and the rapidity of one final-state particle in Fig. 13 The three panels give 
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Figure 14: Distributions for squark pair production pp — > ulUl as a function of the squark transverse momentum (left) and 
rapidity (right). The curves for the central scales we normalize to unity. The scale uncertainty curves we normalize to the same 
central value. The yellow area shows the scale uncertainty, e.g. da/dprifJ- /2) — da / 'dpT(2p°), compared to the purple area 
giving da MLM /dp T - dcr NLO /dp T . We examine the benchmark points CMSSM 10.2.2 and mGMSB 2.1.2. 



LO, NLO, and merged predictions for squark pair production pp — > ulUl- The NLO and the merging histograms are 
normalized to unity, while the LO distribution is shown to scale. As expected, we do not find any kind of significant 
difference between the NLO and the jet merging results nor any correlations between the rapidity and transverse 
momentum. 



B. Scale uncertainties 



The stabilization of the (unphysical) dependence with respect to the choice of the rcnormalization and factorization 
scales is a most prominent feature of NLO calculations. These improvements, which we have already analyzed for the 
total cross section, we reexamine now for the distributions. Again, we study squark pair production pp — ¥ u^u^. In 
Figure [T4| we present the squark transverse momentum and rapidity distributions. 

We first overlay the normalized distributions from the fixed-order NLO calculation (solid, red line) with the central 
scale choice and the jet merging results (dashed, green line) for the CMSSM 10.2.2 and mGMSB 2.1.2 parameter 




Figure 15: K factor as a function of Pt(ml) and y{uL) for squark production pp — > Ulul- The band shows a scale variation 
H°/2 < /i < 2fi°. All MSSM parameters we fix to CMSSM 10.2.2 and mGMSB 2.1.2. 
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points. For the NLO curve we compute the envelope varying the renormalization and factorization scales between 
/i°/2 and 2/i°, keeping the normalization relative to the central scale choice. This should give a realistic estimate of 
the theoretical uncertainty. 

Two differences we show separately: first, the yellow (light) histogram shows the difference da/dpri^ /2) — 
da /dpr{2fx°). It indicates a theoretical uncertainty of 0(10%) on the distribution, with no obvious caveats. In addi- 
tion, we show the difference between the central NLO prediction and MLM multi-jet merging da MhM /dpT—da NLO /dpr 
point-by-point in the purple (dark) histogram. Both comparisons we repeat for the squark rapidity distributions. We 
see that when it comes to normalized distributions the NLO and MLM multi-jet merging predictions are in excellent 
agreement, for example compared to the sizable NLO scale dependence. 

A complementary viewpoint in terms of phase space dependent K factors we display in Fig. |15| The NLO histograms 
using central scales p° are supplemented by a band showing a simultaneous renormalization and factorization scale 
dependence at NLO. We confirm that the K factors remain stable and relatively constant for the transverse momentum 
and the central rapidity regime. From the above discussion we know that the slight change in the K factor over the 
entire phase space should correspond to distributions computed using multi-jet merging. This result we interpret as 
a strong argument in favor of the conventional procedure, where a global iVfactor or event re-weighting to NLO is 
applied to kinematic distributions generated via multi-jet merging. 

C. Non-degenerate squarks 

In Sec. [IT] we discuss the effect of a general squark mass spectrum on the different LHC production rates and find 
that it largely factorizes from the NLO corrections, i.e. the K factors only change at the per-cent level. The impact of 
general squark mass spectra becomes much more apparent at the distribution level. In Fig. [16] we display the squark 
transverse momentum and rapidity distributions. We single out one particular production channel, pp — > u^g and 
examine the following representative situations: (i) mass-degenerate squarks, with m,q = 800 GeV; (ii) a right-left 
splitting Amj;_i = 200 GeV between the right-handed and left-handed squarks; and (iii) a similar down-up splitting 
Arrid-u = 200 GeV. The remaining MSSM parameters we anchor as in the mGMSB 2.1.2 benchmark point defined 
in Table [I] Most importantly, we keep the final-state mass constant, so the differences between these three scenarios 
decouple from the leading phase space effects and instead constitute a genuine NLO reflect. 

The finite mass splitting between squarks induces a shift in the kinematic distributions in the direction of slightly 
harder and more central final-state squarks. We can trace this back to the real emission corrections shown in Fig. |16| 




Figure 16: Normalized transverse momentum (left) and rapidity distributions (right) for squark-gluino production pp — > ul§- 
We assume (i) mass-degenerate squarks with rriq — 800 GeV; (ii) a common mass splitting, Ams_t = 200 GeV; (iii) a common 
mass splitting, Am,j_„ = 200 GeV. The central MSSM parameters we fix as in mGMSB 2.1.2 benchmark. The Feynman 
diagrams to the right describe the squark-gluino fusion mechanism responsible for the significant differences. 
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They describe a fusion mechanism where the bulk contribution arises from internal squark and gluino propagators at 
very small virtuality, i.e. when these particles are almost on-shell. As a result, they become particularly sensitive to 
variations of the squark masses, even if the final-state squarks masses remain unchanged. 

As we can see in Fig. 16 the effect of an 0(20%) mass splitting between up-type and down-type squarks saturates 
the NLO uncertainty on the transverse momentum distributions. Of course, the squark mass spectrum is not a source 
of theory uncertainty which could be captured by the scale dependence. Therefore, it might be useful to estimate its 
effect on LHC analyses independently. 



IV. SUMMARY 



MadGolem is a novel approach to the automated computation of total cross sections and distributions for new 
heavy particles to next-to-leading order. It can be used as an add-on to Madgraph, making use of its interfaces to 
new models as well as to the event generation. In this paper we present a comprehensive overview of supersymmetric 
particle production together with many details of the MadGolem implementation. While MadGolem is not fully 
public yet, a fully functional test version can be obtained from the authors upon request. 

In our application to squark and gluino production we reproduce all relevant Prospino results and extend currently 
available studies in several ways: 

• We evaluate NLO corrections to total rates for a completely general supersymmetric mass spectrum. For 
moderate mass splittings the leading effects of non-degenerate squarks factorize from the LO results, while the 
effect on K factors stays at the level of few per-cent. 

• Instead of identifying the factorization and renormalization scales in the estimate of the theoretical uncertainty 
we vary both scales independently. For heavy strongly interacting new particles the envelope of all possible 
scales agrees with a simultaneous scale variation. 

• Squark and gluino distributions including the full NLO corrections and based on multi-jet merging agree very 
well within the NLO error bands. 

• The effect of non-degenerate spectra on the squark and gluino distributions is clearly visible and can exceed the 
perturbative uncertainty already at moderate mass splittings. 

• The composition of the NLO corrections from different classes of diagrams and with it the dependence of the 
K factors on the mass of the produced particles is significantly different for quark-antiquark vs gluonic initial 
states, i.e. moving from Tevatron to LHC. 



In addition to these specific conclusions we emphasize that with MadGolem this kind of study can be easily 
repeated for any kind of heavy new particles at the LHC. In the light of the available LHC results this might be useful 
not only for general simplified models but also for specific models outside the usual MSSM model and parameter space. 
For this purpose we include an appendix which provides all necessary information on the infrared dipole subtraction 
as well as on a proper on-shell subtraction as implemented in MadGolem . 



Acknowledgments 

We would like to thank Thomas Binoth, who started this project with us and unfortunately cannot see its comple- 
tion. The work presented here has been in part supported by the Concerted Research action "Supersymmetric Models 
and their Signatures at the Large Hadron Collider" and the Research Council of the Vrije Universiteit Brussel and by 
the Belgian Federal Science Policy Office through the Interuniversity Attraction Pole IAP VI/11. DC acknowledges 
support by the International Max Planck Research School for Precision Tests of Fundamental Symmetries. 



17 



Appendix A: Supersymmetric dipoles 



In this appendix we present the unintegrated and integrated dipoles required for SUSY-QCD calculations 26J includ- 
ing a phase space constraint [29]. They are implemented as an independent add-on to the MadDipole package [38] 
and are part of the automated MadGolem framework. 

There exist two major approaches to remove soft and collinear singularities: phase space slicing and phase space 
subtraction [35) . A simple toy example captures their main features and highlights the role of an FKS-like phase 
space constraint [5S]. Let us consider the dimensionally regularized integral J Q dxf(x)/x 1 ~ e with e > 0. Phase space 
slicing based on a small parameter a yields 

rf-l — e I ry.1 — e I rrl — E v J 

q X J a X Jq X 

dx^- + ^- + f(0)loga + O(a;e) . (Al) 
x e 

It requires a choice of small enough a to reach a numerical plateau. A numerically more stable approach is phase 
space subtraction, where the same integral becomes 

Jo x Jo x Jo x 

' fcM e Ma ) + /(0) + oga + 

x e 

Here the divergence is subtracted locally and the final result no longer depends on a. The logarithmic a-dependence is 
exactly cancelled in the total result. This feature we can turn into a sensitive numerical test when varying < a < 1. 



For small values of a we only need to evaluate part of the integrand of Eq.(A2), which speeds up the calculation 



The toy model of Eq.([A2| carries the essence of the Catani-Seymour subtraction method. Real emission of quarks 
and gluons (dcr real ) leads to IR divergences after an integration over the emission phase space. Its regularization 
relies on a local subtraction term (dcr A ) which reflects the universality of the soft and collinear limits. The divergence 
cancels over the same n-particle phase space, 



<5ct nlo = f (dalTo - dat. e=0 ) + / ( d<x virtual + da collincar + f dcrA 



(A3) 



Below, we present the unintegrated dipoles da^ as well as the integrated dipoles da^ including their a dependence. 
They are crucial for SUSY-QCD processes or other NLO QCD predictions beyond the Standard Model. Our extended 
set of massive Catani-Seymour dipoles with explicit a dependence has several practical advantages: 

• tuning a we reduce the subtraction phase space and hence the number of events for which the real-emission 
matrix element and the subtraction fall into different bins; the so-called binning problem. 

• choosing a < 1 we evaluate the subtraction terms only in the phase space region where they matter, i.e. close 
to the IR divergences. 

• our final result should not depend on a. This serves as a test for example of the adequate coverage of all the 
singularities or the relative normalization of the two-particle and three-particle phase space. 

In the MSSM gluino and squark interactions induced by the covariant derivatives glfig, {D^ql 2 give rise to new 
IR divergences which are absent in the Standard Model. The emission of a soft gluon from these particles requires 
new final-final dipoles -Di^fc and final- initial dipoles Dfj. Initial-initial and initial-final configurations can also have 
a squark or gluino as spectator, but the dipole only carries information about the mass of the colored spectator, not 
about its spin. This means we can simply use the massive Standard Model dipoles |26j with an extra SUSY particle 
in the final state. To make this Appendix most useful we will firmly stick to the conventions of Ref. [26) . 



We start with a collection of formulas for final-final dipoles. The expression for the unintegrated dipole is given by 

2ft V T l 
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where |..., ij, k, ...) represents the amplitude for the factorized born process, which in the special case of the SUSY 
dipoles is made by the removal of the gluon from the diagonal splitting q(pij) — >• q(pj)g(pi). The color matrix 
TfcTij/Tfj acts on the born amplitude k, ...) giving the proper color factor. 

To compute the integrated dipoles we integrate over the one-particle phase space [dpi(pij , pk)] with the spin average 
matrices (V^t), according to Eq.(5.22) of Ref. [26] : 

1 a s 1 / Airp, 2 \ c 

[dpi (pij,Pk)\ ~ ; -2 5-(Vi i)fe ) = — — -r — w Iij,k (e) , (A5) 

[Pi+P]) - m fj \ / V / 

where the squark dipole function, (s|V g q k|s') ; i s given by Eq.(C.l) of the same reference, 



8irfi 2e a s CF 



1 - z 3 (1 - yij.k) v ij>k \ PiPj 



5ss< = ^ 3 t k)Ss °' ■ (A6) 



Compared to a massive quark the squark structure is much simpler. This is because for scalars the labels s and s' 
are merely tagging the helicity of the associated quark partners without any effect on the squark splitting. 

The integrated dipole I gq ^ we decompose into an soft or eikonal part 7 clk and a collinear integral 1%-^ evaluated 
in 4 — 2e dimensions, 

I 9 q,k {p q , Pk\ c) = C F [2/ oik fi k ; e) + J^ k {p q , (i k ;e)] 



.,1 / 2 \ 1 9 1 2 

/C1 = 7T l °EP ~ logP log 1 - (M§ + Pk) - - log Pq ~ - log p k 



+ C2 + 2Li 2 (-p) - 2Li 2 (1 - p) - l -U 2 (1 - p|) - ^Li 2 (1 - ^) 

rcoll _ 2 1 2 ..^2 ,.2A , 4/i fc (// fc -l) 



O = ; - ra? - ^ + 6 - 21 °s (a - - 4) + i-7 J ■ (A7) 

e e Pg v ' 1 Pq Pk 

The rescaled masses p, n and the variables p and p„ associated with the splitting ij —> i j and the spectator k are 
defined in terms of the final state momenta pi, Pj and pk as 



Pn 



V(Pi+Pj +Pk) 2 



V / M 1 ^4) 



1 — k . v 
P = 4 / -, , ~ ' Wlth = — j 2 2" 

V 1 + 1 ~ My ~ P\ 

t \ j l-^k + 2pl/(l-p 2 -pl) 
pn{P],Pk) = \ — — ; — , tjz - 2 2T i n = 3, k ), (A8) 

V 1 + + 2/4/ (1 - - Pk) 

with A denoting the Kallen function 

A (x, y, z) = x 2 +y 2 + z 2 - 2xy - 2xz - 2yz . (A9) 
The splitting kinematics we describe using 

Zj = l ■ and y ijk = ■ — >y + = l- g 2 2 ■ ( A1 °) 

PiPk + PjPk P%P 3 + PiPk + PjPk 1 - Pi - Pj ~ P k 

Just like for massive quarks there is no collinear singularity, so the most divergent term in the I gg ,k (e) is a single 1/e 
pole. 

To include the phase space parameter a into the massive squark dipole we limit the dipole function to small values 
of Vij,k/y+ 

n,„,j, ■ /.),,,,,<-> ■ a) ae(0,l]. (All) 
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For the integrated dipole I g q,k (e) we start from Eq.(A7l and subtract the finite term including the same kinematic 
condition as Eq.( All I 



Igq.k (e, a) = I g q :k (e) + Algq.k (a) 

2tt 



gq,k 



(e) 



[ d Pg (PgqiPk)] 



ZPgPq 



e 



Ugq,k 



> a 



(A12) 



The finite part we can evaluate in four dimensions, because by definition there exists no divergence in the region 
Vgq,k/y + >a. The eikonal part 2/[l — z q (1 — y gqt k)\ is the same for (sIVgQ^ls') and (sIV^^s'), so in Eq.(A12| we 
insert Eq.(A7l from our appendix and Eq.(A.9) from Ref. [59"] . 
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The collinear part is different for squarks, so we supplement its form in Eq.(A7l by 

gq 



Ad (a) 



Cf_ 
'2tt 2 



(1 - ii k y - n\ 

2 (1 - a) + log a 



1-/3-/3 



(A13) 

(A14) 

(A15) 
(A16) 

(A17) 

(A18) 



Following the same logic we tackle the final-initial dipoles. The final-initial dipole function is given by Eq.(C.3) of 
Ref. [26], 

(V» 9 ) = 8V £ «^ ^ ry- - 2 - — M ■ (A19) 
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The integrated dipole function I g ~ becomes 

I a gg (x; e) = C F \[3 a g - q (x, N )) + +S(l-x) (j a g f ( W ; e) + J g a f S (/.<?))] + O (e) , 



(A20) 



with the three contributions I g g 



. 1 + log(l - x + uM) 




log (2 + n\ - x) 
log (1 + /if) 2 



4- 



(A21) 



In analogy to the final-final case of Eqs.( All ) and ( A12 1 we introduce a phase space cutoff 



Dgq -> D a g? Q (a - 1 + x gg , a 
0(1- a- x) 



M a q Aa) = -C F 



l-x 



-2 + 2 log 1 



1 + /if - x 



(A22) 



where the kinematic variable Xij. a is given by 



PaPi + PaPj - PiPj + 



111 
m i 3 ~ m i ~ 7Tl j 



% 3-> a 



PaPi + PaPj 



(A23) 



As an example for numerous numerical tests of our dipole implementation, we discuss soft gluon emission off 
the hard process e + e~ — > uru* r . In the left panel of Fig. 17 we show how the dipole subtraction cancels the IR 
divergence locally, i.e. point by point. The numerical agreement of the real emission matrix element with the 
dipole subtraction term improves for softer gluons. In the soft limit both terms grow as l/E 2 g . Even though we find 
l-^rcai ~ E D g q t k\ ~ 1/E g the phase space factor E g dE g cancels this dependence. 

dipolcs 

In the right panel of Fig. [17] we show the a dependence for the final-final squark dipole. Both, the real emission 
and the integrated dipole, depend separately on a. Their sum is numerically stable over many orders of magnitude 
and down to a = 10~ 9 . 



Appendix B: On-shell subtraction 



From single top production it is well known that when including NLO corrections we have to avoid double counting 
of diagrams which are attributed to different physics processes. As an example we consider real emission corrections 
to squark pair production pp — > qq: the partonic sub-channels with an additional quark in the final state qg — ¥ qqq 
display a peculiar behavior which we illustrate in Fig. 18 The diagrams (a) and (b) are part of the genuine NLO 



corrections to squark pair production. In contrast, the diagrams (c) and (d) can be interpreted in two ways: 



qg -> qg x BR(g -> qq) 



squark pair production 
squark-gluino production 



(Bl) 



The first interpretation simply assumes NLO corrections to the hard process pp(qg) qq and is generally valid 
for on-shell and off-shell gluinos. The second interpretation points to the LO process for qg — > qg followed by the 
branching BR(g — > qq) and implicitly assumes an on-shell gluino. For a mass hierarchy m g > m g we can therefore 
separate the two assignments into off-shell and on-shell gluinos. This distinction avoids double counting and is the 
basis of our on-shell subtraction scheme. Approaches to tackle this problem include 
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Figure 17: Left: real emission matrix element (red circles) and the dipole subtraction (black crosses inside) towards the soft 
limit y 9 q t k — > 0. Right: a dependence for final-final squark dipoles. 



• a slicing procedure which separates the phase space related to the on-shell emissions and removes the on-shell 
divergence by requiring \y/Sqq — m§| > 8 [ID]. Phase space methods of this kind do not offer a cancellation of 
the 8 dependence and do not act locally in phase space. Moreover, as a pure phase space approach it does not 
allow for a proper separation into finite, on-shell and interference contributions which is crucial for a reliable 
rate prediction. 

• diagram removal where the resonant diagrams are removed by hand. In lucky cases this method might work 
in the limit T/m <C 1 |41j . but it ignores any kind of interference contributions which do not actually have to 
vanish in narrow width limit. This scheme is theoretically poorly motivated in many ways. 

• local on-shell subtraction in the so-called Prospino scheme [TTJ [27] which under the name 'diagram subtraction' 
is also used in the single-top computation of Mc@nlo. This method used in MadGolem . 

To define the on-shell subtraction we split the contributions of the matrix element in two parts: the first piece 
concerns the resonant diagrams (c) and (d) and is denoted as A^res, while the second piece represents the non- 
resonant (remnant) diagrams (a) and (b) as .M rem . Note that this separation is defined at the amplitude level and 
not based on the amplitude squared. The full matrix element squared becomes 

\M\ 2 = \M rcs \ 2 + 2Re(M* es M rem ) + \M Icm \ 2 . (B2) 




Figure 18: Sample diagrams for the real-emission corrections to squark pair production with an additional quark in the final 
state. 
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Figure 19: Kinematic variables for the on-shell subtraction. 



The divergent propagator in A^ res we regularize as a Breit-Wigner propagator 

1 1 



Pii ~ ™% 



(B3) 



where my is the mass of the mother particle in the splitting ij — > i j, as shown in Fig. 19 and is a regulator. 



As explained above, a possible double counting is limited to the on-shell configuration in |.M r0 s| 2 . To remove it we 
define a local subtraction term da and include it in complete analogy to the Catani-Seymour dipole subtraction 
Eq. ( A3 1 , such that the total cross section is given by 



5a 



NLO 



n+1 



(dal°S-da^ e=0 -da™ ) + 



da" 



rtual 



da' 



collii 



dai 

1 / e=0 



(B4) 



The extra subtraction term da os is |A4 res | 2 with its momenta remapped to the on-shell kinematics, 

1 



da os — Q(S — (rriij + rrik) 2 ) 0(m.y — rrii — nij) 



(p 



1 



m 2 T 2 - 

i 3 



(B5) 



remapped 



The kinematic configuration is depicted in Fig. 19 The two step functions in Eq.(B5) ensure that the partonic 



center-of-mass energy is sufficient to produce the intermediate on-shell particle and that it can decay on-shell into the 
two final state particles. The ratio of the Breit-Wigner functions ensures that the subtraction has the same profile as 
the original |7W res | 2 over the entire phase space. In the small width limit this ratio reproduces the delta distribution 
which factorizes the 2 — » 3 diagrams into a x BR. 

Note that this method works with a mathematical regulator Fy which can be related to the physical width as 
in the Mc@nlo implementation; alternatively we can go into the well-defined limit Fy <C rriij used in the original 
Prospino implementation. 

In summary, this on-shell subtraction implemented in MadGolem exhibits several attractive features when it comes 
to prediction of total and differential cross sections. First, it subtracts all on-shell divergences point-by-point over the 
entire phase space. This means that all distributions are automatically safe. Second, it preserves gauge invariance at 
least in the narrow-width limit. Third, it takes into account spin correlations, because it includes the full 2 — > 3 matrix 
element. Fourth, it keeps track of the interference of the resonant and non-resonant terms, 2Re(A^* cs A4 rC m), which 
can be numerically sizable. Finally, Fig. 
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and the narrow- width limit Yij/m.ij — > 



shows that it smoothly interpolates between a finite width Tij /rriij ~ 0.1 



Appendix C: One-loop amplitudes 

The virtual corrections module of MadGolem is based on the Feynman-diagrammatic approach defined in Ref |24j . 
It calculates the virtual corrections to any 2—^2 process such that it can be applied to the Standard Model, MSSM, 
or any other renormalizable theory. 

The module uses Qgraf [23] , FORM j52J, Maple, and the Golem95 integral library [2S] to produce Fortran90 
code that calculates the virtual cross section for a 2 — > 2 process, given a set of phase space points. It also produces 
analytical Maple output in the form of partial amplitudes. Each virtual diagram is broken down according to its 
color structure, helicity, and scalar integrals. This allows for a careful test before the numerical calculation is even 
started. The approach can be divided into four main steps: 
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Figure 20: NLO contributions from intermediate on-shell particles in uu — > UrXi + X production as a function of Tn/m-a. The 
squark width acts as a cutoff in the Prospino subtraction scheme [111 I27| . The masses are chosen to illustrate all different 
resonant channels; virtual corrections are not included. 



1. generate the tree level diagrams, counterterms, and one-loop diagrams with Qgraf and translate the output into 
Form code suitable for symbolic manipulation. The analytic structures keep track of external wave functions, 
vertex couplings and internal propagators, color factors, Lorentz structure, and the overall sign from external 
fcrmions. 

2. map the analytic structures onto partial amplitudes using a basis in color, helicity and tensor structures using 
the spinor-helicity formalism. 

3. apply an analytical reduction of tensor integrals to scalar loop integrals based on the Golem reduction 
scheme [24 1 125 ] , 

4. collect all results and insert the renormalization constants into the counterterms. The final output is available 
as analytical partial amplitudes in Maple and as numerical Fortran90 output. The latter is implemented 
into the Madgraph structure. 

These four blocks are coordinated and run by the Perl script run_golem.pl. We will describe them in detail below. 

In the first step Qgraf generates all Born diagrams, counterterms, and QCD one-loop corrections for a hard process 
specified in a MADGRAPH-like file proc_card.dat. The topological rules Qgraf obtains from modified MadGraph 
model files. In addition to the familiar MadGraph options, we include novel functionalities specific to a NLO 
calculation; for instance, the flag nlotype in the process card generates either pure QCD (gluon mediated), or full 
SUSY-QCD virtual corrections. This division relies on a constrained set of propagators within Qgraf. Note that 
SUSY-QCD corrections also include loop diagrams which do not involve either gluons or gluinos, e.g. mediated by 
the four-squark vertex. Therefore, Qgraf first includes all loop diagrams with gluons, gluinos and squarks as well as 
Faddeev-Popov ghosts. The number of loop diagrams is reduced later by checking the order of a s . 

Counterterms are generated automatically by Qgraf via tree-level diagrams containing place-holders for all renor- 
malization constants. These renormalization constants depend on 0(a s ) corrections to a set of two-point functions 
involving the different colored particles present within a given model. We provide them as a set of separated (model- 
dependent) libraries, implemented as Maple list files. 

Additional topological constraints, e.g. requiring only gluonic i-channel contributions, or only self-energy correc- 
tions, can be included via ruri-golem.pl At this stage the code groups topologically equivalent structures and applies 
loop filtering techniques e.g. removing diagrams which are trivially zero. The Perl script also assigns the overall 
sign to each diagram, because standard Qgraf is not valid for Majorana fermions. 

Once the Qgraf output is filtered by run_golem.pl, the remaining set of diagrams is processed in Form to apply 
Feynman rules. Assigning the correct fermion flow is crucial for diagrams with Majorana particles |43j . 

The second step of MadGolem treats the color and helicity structure of the Qgraf and Form output. For the 
QCD structure of each Feynman diagram FORM uses a color-flow decomposition [33]. Each external gluon is matched 
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with an adjoint generator T?-, which means we can re-write the gauge structure using delta functions in color-space, 
and factorize it from the remaining amplitude. This way the color flow within the amplitude becomes more apparent. 

The spin structure of each diagram is also manipulated in FORM. Using the spinor-helicity formalism [151 ES] the 
amplitude is projected onto a set of helicity amplitudes. Each fermion pair and vector boson is re-written as massless 
spinor products, of which we take the traces. This way each diagram is expressed in the Mandelstam variables s, t 
and u, with a spinor product prefactor. Massive spinors require a helicity projection in the direction of an auxiliary 
reference vector which we choose to be the light-like momentum of one of the other external particles fc,. Whenever 
this is not allowed we instead use k$ = (E% + i?2)(l, 1,0, 0)/2. In this case the kinematic structure is no longer 
defined by the usual Mandelstam variables, so we also include S15, S25, S35, S45, and ie^pa-k^k^k^k^ . This additional 
reference vector seriously impacts the tensor reduction described below and slows down the amplitude generation and 
evaluation. 

In step three of MadGolem we simplify the loop diagrams using the Golem approach [35] ■ For a fully 
analytical reduction of tensor integrals to a linear combination of scalar integrals we rely on a combination of Form 
and Maple. All one-loop integrals are regularized dimensionally; internal momenta and gamma matrices are split 
into four-dimensional and (— 2e)-dimensional components, with the latter only contributing at 0(e) [47] . Tensor loop 
integrals we simplified using a Passarino-Veltman reduction [48j . This reduces an TV-point tensor integral of rank r to 
a scalar A-point integral plus a series of integrals with fewer external legs and reduced rank. The final result we can 
express in terms of known scalar integrals (-Do, Co, Bq, Aq) in 4 — 2e dimensions. Their divergence structure is simple: 
IR poles arise purely from D and Co, UV poles arise purely from B and Aq. The only exception to this rule is the 
scalar two-point function 

W-W?^ ^7^ = 1^?^ (-"-)' < C1 » 

J (2it) u q z (q + p) 1 16tH F(l - e) V e UV eiR/ 

which in this schematic notation is IR and UV divergent. Its finite part vanishes, but in our calculation we need 
to keep track of its IR and UV poles separately. As mentioned above, four-point tensor integrals with k§ ■ I in the 
numerator (I standing for the internal loop momentum) cannot be simplified using the Passarino-Veltman approach 
and are kept as un-reduced form factors to be numerically processed by Golem95. 

In the final step we collect all partial amplitudes for a given process using Maple. Two analytical output files 
contain all information about the Born amplitude and the renormalized virtual amplitude: 

• AMP_TREE.mapout lists the total non-zero Born amplitudes, sorted by diagram, helicity, and color represen- 
tation. If the flag nlosymsimp is enabled in the MadGraph process card, the helicity amplitudes are tested for 
the possible symmetry = A4^ x I*, where {A'} is a different helicity from {A}. Only the minimal set of 
helicity amplitudes is kept, along with a note of which helicities are conjugates of which. This greatly reduces 
the size of the output for pure QCD or QED processes. 

• AMP_LOOP.mapout lists all finite loop amplitudes as kinematic coefficients sorted by diagram, helicity, color 
representation, and type (scalar integral, form factor, or number). In the same format it also lists the countert- 
erms after the renormalization constants have been inserted. The simplification flag nlosymsimp is also applied 
to the loop amplitudes. 

For a numerical evaluation we do not rely on this analytic output. Instead, MadGolem writes several Fortran90 
routines for the computation of the virtual corrections: 

• libcoeffs-all.so and Ubcoeffs_alLtree.so contain the amplitude coefficients for the virtual corrections and Born 
amplitudes. For size reasons we generate a separate library for each partial amplitude. Each library we pre- 
compile before linking them dynamically and launching them at runtime. 



golem(k,mu, amplitude-array) takes the external four- momenta &i ,2,3,4 and the renormalization scale p, and re- 
turns 



amplitude -array = ^a , — -, — y-, |A4|| orn , auv) • (C2) 

The different aj are defined through the interference term between Born and virtual amplitude, 

|M|i loop = 2Re (M* hoIn M viTt ) = a + — + ^ , (C3) 

and correspond to the finite contribution (ao), and the coefficients of the single (a_i) and double (0-2) IR poles. 
The Born term is included for comparison with MadGraph, and ayy returns the numerical value of the UV 
pole which is zero after proper renormalization. All results are averaged over initial state colors and helicities. 
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• virtuaLcorrections.f90 contains the subroutine golem(k,mu,amplitude_array), which calls the integral library 
Golem95 and the coefficient libraries libcoeffs*.so. It calculates the fully renormalized matrix element at one- 
loop level. The three values Uq, a— i, 0,-2 a- re combined with the integrated dipoles from the real emission 
corrections for the complete NLO corrections to a 2 — > 2 process. The cancellation of the single and double 
poles is automatically checked in our numerical implementation. 

As is appropriate for an automized NLO tool we have undertaken an exhaustive program of checks to ensure the 
robustness and reliability of our MadGolem . We have calculated the (SUSY)-QCD one-loop corrections for a large 
set of 2 — > 2 processes both in the Standard Model and the MSSM, covering all representative possibilities of spins, 
color charges, interaction patterns and topologies. The cancellation of the UV, IR and OS divergences, as well as 
the gauge invariance of the overall result, can be confirmed numerically, for some specific cases also analytically. 
The finite renormalized one-loop amplitudes we have systematically compared with FeynArts, FormCalc and 
LoopTools [35] . 



Appendix D: Renormalization 



As discussed in the previous appendix, we automatically generate the ultraviolet counter terms using the tree-level 
amplitude from Qgraf. As an input we express all renormalization constants in terms of two-point functions as a 
separate library. The current MadGolem implementation fully supports renormalized QCD effects for the Standard 
Model, the MSSM, sgluons [T7], and other new physics models. To document our notation we give all relevant 
expressions for the renormalization of supersymmetric QCD here. 

The renormalization constants we define through the relation between the bare and the renormalized fields, masses 
and the coupling constant: 



(0) 



4 /2 * 



(0) 



mq, + 5 my 



9. 



(0) 



9s + 5g s 



(Dl) 



The field \fr = q,q,g,g denotes all strongly interacting MSSM fields. We express the SUSY-QCD counter terms to 
vertices and propagators in Table [TV] 



The actual counter terms, presented below, we include in a separate library. The strong coupling constant we renor- 
malize in the MS scheme and explicitly decouple all particles heavier than the bottom quark. This zero-momentum 
subtraction scheme I50[ 15 lj leaves us with the renormalization group running of a s to light colored particles only. 
It corresponds to the measured value of the strong coupling, for example in a combined fit with the parton densities. 
Its renormalization constant reads 



<5 g s 
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The UV divergence appears as 1/e = (47r) e /r(l — e) = 1/e — Je + log(47r) 
colored particles contribute to the beta function 



0(e). Both light (L) and heavy (H) 



A) = /?o L + Po 
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C A 



(D3) 



MadGolem sets the number of active flavors to nf = 5. The SUc(3) color factors are Cf = 4/3 and Ca = 3. 



The analytic form of all renormalization constants we reduce down to one-point and two-point scalar one-loop 
functions, which we handle by means of the standard 't Hooft-Vcltman dimensional regularization scheme in 4 — 2e 
dimensions. The field and mass renormalization constants we compute from the one-loop self-energies which involve 
virtual gluons and gluinos. All fields are renormalized on-shell. In addition, for the gluon field we subtract the 
heavy modes consistently with our g s renormalization scheme. The underlying Slavnov- Taylor identities link the 
corresponding finite counter terms as 5 Zq = —25g s . 

In addition, we need to pay attention to dimensional regularization breaking supersymmetry through a mismatch 
of two gaugino and the 2 — 2e gauge vector degrees of freedom [3D]. As a result, the Yukawa coupling g s appearing 
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iS Zg — nig 8 Zg — 5 VClg 



Table IV: Strong interaction counter terms for the MSSM. The finite supersymmetry-restoring counter term <5susy is given in 
Eq.fD4l. 



in the qqg vertex departs from g s . We restore supersymmetry by hand, forcing g s — g s . The corresponding hnite 
counter term can be computed using dimensional reduction, 



9s a s ( 2 3 



6 S 
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3 4tt 



(D4) 



Finally, we quote the analytical expressions for the field and mass renormalization. For the scalar one-point and 
two-point functions we adopt the notation of Ref. [52] . The corrections to the massless quarks including the non-chiral 
SUSY contributions are 



Bo(0,0 J Q)+B o (0,m| > m| ) 



(D5) 



The corresponding squark fields {q s =L/ii) an d mass are renormalized as 



Z7T 
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The gluon wave function renormalization, linked to the counter term for the strong coupling, is 



77 ^g 
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Finally, the gluino field and mass renormalization constants are 
a s 



SZ- n 



■Ca 



l + 4m 2 9 ^(m|,0,m|)- Ao(ro ^ 



+ 



47T 

E [A (m?)-(m?+m?)S (ml,0 t m?)-2m?(m|-m|)^(m| ) 0,m?)] 

^ light (s)quarks 

a s To 2/2 2 2\ r>l / 2 2 2 \ , / 

2^ [ 2 m s {m s -m q -m s )B {m s ,m q ,m Sa ) + {-t 

9 heavy (s)quarks 



2 2\ n / 2 2 2\ 

m s - m 5 ) B (m d ,m q ,m s ) 



Oh 

n 



6 m- 



47T 

Q-5 



heavy (s)quarks 

C A m„ 



+ A Q (m\)- A (m 2 q )] 
iu;, m„ R q sl R q s2 B'Q{m g ,m 2 q ,m qs ) 



1 + 3- 



87T7T1S 



J] [A Q (m?) + (m| - m?) B (m?, 0, m§)] 



light (s)quarks 



8lT Win 



53 [ A o(m 2 q) - A (m 2 q ) - (m\ -m\- m~) B (m 2 ,,m 2 q ,m 2 q )] 



Q 

heavy (s)quarks 

^ 5Z m qRllRl2 B o{ m li m % m l s )- 

heavy (s)quarks 



(D7) 



(D8) 



The sum over heavy squarks covers all squark flavors corresponding to heavy quarks. We usually consider the bottom 
quark massless, which means that only the two stop eigenstates feel top mass effects. However, the bottom/sbottom 
loops can be trivially moved from the light to the heavy category. The stop mass eigenstates tip are related to the 
electroweak interaction bases through a rotation with R = ±1. 
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